Multifractal properties of power-law time sequences; application to ricepiles 
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We study the properties of time sequences extracted from a self-organized critical system, within 
the framework of the mathematical multifractal analysis. To this end, we propose a fixed-mass al- 
gorithm, well suited to deal with highly inhomogeneous one dimensional multifractal measures. We 
find that the fixed mass (dual) spectrum of generalized dimensions depends on both the system size 
L and the length N of the sequence considered, being however stable when these two parameters 
are kept fixed. A finite-size scaling relation is proposed, allowing us to define a renormalized spec- 
trum, independent of size effects. We interpret our results as an evidence of extremely long-range 
correlations induced in the sequence by the criticality of the system 
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I. INTRODUCTION 

Self-organized criticality (SOC) has been the subject 
of a great deal of interest, since its introduction by Bak, 
Tang, and Wiesenfeld Q. The main feature of SOC sys- 
tems is that they evolve, driven by means of an external 
force, into a critical state characterized by the absence 
of any characteristic time or length scale. The result- 
ing extremely long range correlations show up through 
the peculiar "1//" power spectrum and the geometrical 
fractal structure. SOC behavior has been observed in 
many cellular automata models of sandpiles Q| , invasion 
percolation 0, biological evolution depinning in ran- 
dom media [Q, and also in some natural systems, such 
as earthquakes jf|. Even though the first cellular au- 
tomaton displaying SOC was conceived to represent the 
dynamics of a sandpile Q , the experimental evidence in- 
dicates that this is not actually the case: Real sandpiles 
are not in a self-organized critical state |]-||- Recently, 
however, Frette et al. || reported SOC behavior in a real 
granular system, a one dimensional ricepile. For grains 
of rice with a considerable aspect ratio, the pile behaves 
critically, this fact being accounted for by the increased 
friction, which is able to counterbalance the inertia effects 
predominant in real sandpiles. 

In a subsequent paper, Christensen et al. [ [l0[ analyzed 
the transport properties of individual grains inside a sta- 
tionary ricepile. They measured the transit time of in- 
dividually colored grains of rice (tracers), defined as the 
time necessary for a grain to escape from the pile. Chris- 
tensen et al. found that the distribution of transit times 
follows a truncated power law form, and that the average 
transport velocity of the grains diminishes as the system 
size increases. A cellular automaton model of a ricepile 
was proposed (the so-called Oslo model) 10 llj], repro- 
ducing the phenomenological behavior of the actual ex- 
periments. Boguna and Corral jl2j have also suggested 
a theoretical scenario for the Oslo ricepile, based on a 
continuous time random walk model. 

The main results of the Oslo experiments and sim- 



ulations can be expressed through a single function, 
the probability distribution of transit times P(T, L)dT, 
which is defined as the probability of a given tracer 
spending a time between T and T + dT inside a pile of 
linear size L. It was found that 



P(T,L)~L- U , T<V, 
P(T,L) ~T~ X , T>L V , 



(1.1) 



where v and x are certain characteristic exponents. The 
experiments provided the values v = 1.50 ±0.20 and x — 
2.40 ± 0.20, whereas the cellular automaton model ren- 
dered the exponents v = 1.30 ±0.10 and x = 2.22 ±0.10 
jn^jOJ . This numerical evidence can be summarized in 
the finite-size scaling ansatz 



P(T,L) = L~ l3 f 



T 
TP 



with 
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\x-x 



for x < 1 
for x > 1 



(1.2) 



(1.3) 



Given that x > 2, and provided that the probability 
distribution is normalized, we have that (3 ~ v and the 
average value of T is finite, (T) ~ L v (see Appendix). 
The fact that x < 3 implies, however, that the second 
moment of the distributio n is infinite, (T 2 ) = oo. 

The finite-size scaling (1.2) compacts the experimen- 
tal data into a useful relationship, which on its turn al- 
lows one to extract valuable conclusions about the sys- 
tem. However, it is actually quite obvious that it is pos- 
sible to extract more information about the ricepile from 
the sequence of transit times, apart from its distribution 
function. In order to gain a different insight into the 
problem, we propose to consider the output of the exper- 
iment from a different point of view. Let us define the set 
S(N, L) as follows: Throw a tracer grain in a stationary 
pile ]l3j of linear size L and measure the time elapsed un- 
til it emerges outside within any avalanche. Performing 
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the same measurement for N different grains, consecu- 
tively thrown in the pile, we can construct the sequence 
S(N,L) = {T n } n=1 pf, where T n is the time, mea- 
sured in units of added grains (the slow time scale fL0[), 
spent inside the pile by the n-th grain, in the sequence 
of N consecutive throws. The set S(N,L) can be in- 
terpreted as a discrete time sequence, assigning to the 
instant n = 1,...,N the value T n . In the Fig. 1 we 
have represented such a sequence, for the transit times 
recorded in the cellular automaton model of ricepile de- 
scribed in Ref. jl(|. The system size is L = 100. Fig. 1(a) 
shows a record of 90000 transit times, whereas Fig. 1(b) 
depicts the 5000 points closer to the center of the pre- 
vious graph. These plots show rather conclusively that 
not only the distribution of transit times is scale-free, but 
also that their sequence is in some sense self-similar. 
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sideration. When computing the multifractal spectrum 
of the sequence S(N,L), we observe that it shows con- 
siderable size effects: The spectrum of generalized di- 
mensions D(q) (to be defined later on) depends on the 
system size L and, even worse, on the sequence length 
N. This fact seems to doom any effort to describe a sin- 
gle well-defined spectrum. However, by analyzing D{q) 
in the limit q — > oo, we observe a power law dependence 
on N and L. Extending the scaling to the whole range 
of q allows us to define a "renormalized" spectrum, truly 
independent of size effects. We interpret our results as an 
effect of the extremely long ranged correlations present 
in the sequence, correlations induced by the criticality of 
the ricepile. 

We have organized this manuscript as follows: In 
Sec. II we review the multifractal analysis of general 
mathematical measures, stressing the difference between 
fixed-size and fixed-mass formalisms. In Sec. Ill we de- 
velop in particular the formalism needed to deal with a 
discrete time sequence. Sec. IV analyses different syn- 
thetic uncorrelated random time sequences. First we 
check the accuracy of the algorithm against sequences of 
known spectra. Then we study a power-law distributed 
random signal, mimicking the real transit time sequences. 
Section V deals with our final goal, actual sequences 
of transit times from numerical simulations of the Oslo 
model. Finally, our conclusions are discussed in Sec. VI. 



II. MULTIFRACTAL ANALYSIS: FIXED-SIZE VS. 
FIXED-MASS FORMALISM 



12 



10 - 



H 
o 



6 - 



50000 



b) 



51000 



52000 53000 

n 



54000 



55000 



Fig. 1: a) Sequence of transit times for 90000 tracer grains 
in a computer simulation of the Oslo ricepile; system size 
L = 100. b) Zoom close to the central section of the previous 
picture. 
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Loosely speaking, we call multifractals 
mathematical sets which can be decomposed into an infi- 
nite set of interwoven subfractals, labelled with an index 
a, each of them characterized by a different fractal di- 
mension /. The collection of these dimensions form the 
so-called multifractal spectrum f(a) [ p^[ . Strictly speak- 
ing, however, it is only possible to assign mathematically 
meaningful multifractal properties to a measure (mathe- 
matical or physical) defined over a given support Jl9| . A 
multifractal measure is completely specified either by its 
multifractal spectrum /(a) or by its spectrum of gener- 
alized dimensions D(q). 

In this section we review the main mathematical defi- 
nitions and properties of multifractal analysis. 



A. General definitions 



In this paper we will extract more information from 
the Oslo ricepile model studying the sequence of transit 
times S(N,L). The method we have employed is that 
of multifractal analysis (which, on the other hand, is not 
new in the field of SOC HQ). To this end, we have 
developed an algorithm particularly well suited to deal 
with one dimensional measures, like the ones under con- 



Following Ref. |T^] (see also p0|), consider a normal- 
ized measure \i defined on a support K C IR d , fJ.(K) = 1. 
Let A be an arbitrary partition of K in non-intersecting 
elements A,, that is, 



K C (J A, and A 4 Q A^ 



i?3, (2-1) 



2 



and let pi and £j, i = 1, . . . , N be the variables that rep- 
resent the weight factor and the size factor corresponding 
to the element A,-, respectively. We define the function 



JV 



(2.2) 



where g and r are any real numbers. The sum runs over 
all the N disconnected parts in which we decompose the 
support of the measure, and the brackets stand for an 
average over different realizations of the measure. For 
any measure, either deterministic or experimental (non 
deterministic), we will assume that, for fine enough par- 
titions, the function < &a(9, t) collapses onto a single con- 
stant value, that is, 




const. 



(2.3) 




Expression ( |2.3[ ) is an implicit equation, allowing one 
to determine r(q) for a given q, or conversely, q(r) 
for a given r. If we assume a partition A in which 
£i = e = const., then the size factor e can be factorizcd 
from the former expression, yielding 



(2.4) 



where N(e) is the number of parts of size s, contain- 
ing a certain measure pi different from zero. From 
this last expression we can compute the function r(q) 
and the generalized dimensions D(q) |f2l"| , ^2f , defined by 
D(q) = r(q) / (q — 1) . The f(a) spectrum is given by the 
Legendre transformation f(a) = imii q {qa — (q— l)D(q)} 
p9| , p3[ . This approach corresponds to the so-called fixed- 
size multifractal formalism (FSF). 

On the other hand, we can select a partition A in which 
Pi = p — const., which yields to 



(2.5) 



where N(p) is the number of parts of measure p, with a 
certain size £j different from zero. From this expression 
we can calculate the function q(r), and then, inverting 
it, compute the spectrum D(q). This second approach 
corresponds to the so-called fixed-mass multifractal for- 
malism (FMF). 

Both FSF and FMF are completely equivalent. In 
order to stress this correspondence, we define the new 
parameters q* = — r and r* = —q and substitute into 
(|2~5l). Now both equations (EO) and (EE" 



Eq. fl2.5|). Now both equations fl2.4| ) and ( |2.5D read the 
same, the only difference being the change of role of pi 
and £i. The equivalence between both formalisms is ex- 
plicitly illustrated by the identities 



D*(q*) = 



q* = -(q-l)D(q), 

q 



(2.6) 



\ + {q-l)D{qy 
with D*{q*) = T*{q*)/{q* - 1). 

B. Box-counting algorithms 

The most common operative numerical implementa- 
tions of multifractal analysis are the so-called fixed-size 
box- counting algorithms |18| . For a given measure n with 
support K C H , they consider the partition sum 



Ze(q) = 



E 



KB) 



(2.7) 



q S 1R, where the sum runs over all the different non- 
empty boxes B of a given side e in a e-grid covering the 
support K, that is 



(2.8) 



fe=i 



Ik being integer numbers. The generalized fractal dimen- 
sions of the measure are defined by the limit 



W q-le^O loge 



(2.9) 



and numerically estimated through a linear regression of 



1 



log Z £ (q) 



(2.10) 



against loge. 

Wi thin this formalism, the mathematical defini- 
tion ( 2^9 ) is strictly valid for positive q J24J]. Numerical 



estimates work well for q > 1 in d > 2, and render usually 
incorrect results for q < J25|-|27|]. This fact is obviously 
due to the presence of boxes B with an unnaturally small 
measure, which contribute to the function Z with diverg- 
ing terms. In those cases, one is forced to apply different 
prescriptions [p7 28 1. 

The box-counting version of the fixed-mass formalism 
is in general harder to implement in d > 1 spatial di- 
mensions. The difficulties reside in the proper selection 
of boxes with a given fixed measure. (For an applica- 
tion in d = 2 see Rcf. f20|| .) From a numerical point of 
view, it is well known that the FMF is a good estimator 
of ge ner alized dimensions for q < (that is, q* > 0, see 
Eq. (pTej) ) and bad for q > (q* < 0). The explanation 
of this behavior is related to the space distribution of the 
measure. The FSF operates well in the dense regions of 
the support, whereas the FMF is specially appropriate 
to deal with its sparse regions. 

As we will see in the next section, however, a fixed- 
mass algorithm is particularly simple to implement for 
one dimensional measures, such as time sequences. 
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III. MULTIFRACTAL FORMALISM FOR 
DISCRETE TIME SEQUENCES 

Fractal geometry and multifractal analysis are well- 
known tools for the study of complex time signals (see for 
instance f2!||3(J and references therein). In this section 
we will specialize the box-counting multifractal analysis 
sketched above for the particular case of a discrete one 
dimensional time sequence. 

We define a general discrete time sequence T(N) as 
any set of N positive real numbers, T(N) = {t n } n= i ni 
t n € 1R + . At this level we will not make any assumption 
about the possible correlations of the sequence. However, 
we will assume that it is the outcome of some physical 
process in a stationary state, and that we can obtain se- 
quences as long as it might be required. 



A. Fixed-size algorithm 

In order to study the multifractal properties of a se- 
quence T(N), we must first provide a meaningful physi- 
cal measure on it. As a first ansatz, we define the naive 
measure fi on the support ]0, N] C H over which the se- 
quence is defined. This measure assigns to a given box 
in ]0, TV] a weight proportional to the sum of the value t n 
of all the points n inside the box Namely, if B(x, e) 
is a ball with center in x and diameter e, then 



Q(N) ^ 



in • 



(3-1) 



where Q(N) = X^=i is a normalization factor such 
that n(]0, N}) = 1. In order to compute the generalized 
dimensions D(q) of fx, consider a partition of ]0, N] into 
boxes of diameter r, in a number N/r, defined by 



B k . r =}{k-l)r,kr], fc = l, 
The partition sum will then read 



,N/r 



(3.2) 



N/r N/r 

^) = EW= P E( E 

k=l ^ y ' k=l (k-l)r<n<kr 



(3.3 i 



The generalized dimensions arc defined through 

l°g | Q(N)i Sfc=l ( J2{k-l)r<n<kr } 



D(q) 



1 



lim 



q — 1 r/JV->0 



log^ 



(3.4) 



The role of e is now played by the reduced diameter of 
the boxes r/N. Numerically, we will obtain an estimate 
of D(q) as the slope of a linear regression of 



N/r 



fc=l (k-l)r<n<kr 



(3.5) 



against log(r/JV). Note that we have dropped the nor- 
malization factor Q(N) q , since it does not depend on r 
and therefore plays no role in the regression. Moreover, 
the elimination of this factor results in general in a better 
performance of the numerical algorithm, except for those 
values of q very close to 1. 



B. Fixed- mass algorithm 

In order to define a fixed-mass algorithm for a discrete 
sequence T(N) , we start by constructing an approximate 
Cantor set Cr{N), composed by a collection of N dis- 
crete points on the interval ]0, 1]. We define the dual 
measure fi* by associating a mass distribution to this 
approximate Cantor set. The distribution corresponds 
to just assigning a mass unity to each one of its points. 
Consider thus the sequence T(N) — {t n } n=1 N , with 

Q(N) — J2n=i an d l et us define the Cantor set Ct{N) 
by 



C r {N) = {x n | < x n < l,n=l, 



with 



1 " 

ivtE^- 



Q(N) 



,N} (3.6) 



(3.7) 



k=l 



We define a measure on Ct(N) through the density func- 
tion 



Pel 



1 N 



(3.8) 



where 6 is the usual Dirac delta function. The measure of 
a ball of center x and diameter e, B(x, e) = ]x — |, x+ |] 
is given by the integral 



(i*(B(x,e)) 



p c (x)dx 



(3.9) 



and is equal to the number of points from Ct{N) con- 
tained in the interval B{x,e). It is easy to verify that 
the dual measure fj,* has holes of finite size: Consider a 
given t p and 



£ < 



Q(N)' 



and define 



X — Xp — i ~r 



2Q(N) 



(3.10) 



(3.11) 
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Then n*(B(x, e)) — 0. If t p is very large, then it will cor- 
respond to a large hole in Ct(N), with a diameter gran ' 
This implies that the fractal dimension of the support 
of 11* would be in general lesser than 1. These regions 
of zero dual measure are related to the regions of large 
naive measure. 

We define the FSF multifractal spectrum of fi*, 
D*(q*), through the partition function Z*(q*) 7 which on 
its turn is defined onto the basis of a set of disjoint in- 
tervals covering ]0, 1]: 



B he =](k-l)e,ke], 



(3.12) 



that is, 



l/e 



k=l 



l/e 



— E 



fc=i 



k-z 



p c (x)dx 



(fc-l)e 



(3.13) 



The generalized dimensions are mathematically defined 
by the limit 



D*(q*) 



i 



q* — 1 £^o loge 



(3.14) 



and numerically evaluated as the slope of a linear fit of 
1 



q* - 1 



log W) 



(3.15) 



against log e. We will drop again the normalization factor 
2V«\ ' 

From a mathematical point of view, this construction 
represents a practical implementation of the notion of in- 
verse multifractal measure discussed by Mandelbrot and 
Riedi in Ref. [flj. Let us show that ji* indeed corre- 
sponds to the inverse of the naive measure defined on the 
original sequence. Consider a box of size Sk, which 
contains points from the Ct(N), and therefore has an 
associated dual measure 



nk_ 
N' 



(3.16) 



Consider that those rik points are the consecutive 
points xi,xi + i, . . . ,xi +nk _i. Assuming that the extreme 
points coincide with the extremes of the interval, then 
we have x;+ nfc -i — xi ~ e k . If we recover the former 
definition of x„ , then 



j i+nfe— 1 ^ / ^ i+nfc— 1 

^^-'-^ot E tj "ow? t, = w ? is "^ fc) ' (3 - 17) 

S — 1 S — 1 S — f — ] — 1 

where Bk is a certain box, associated with the naive measure, with diameter ik ~ Uk/N. Then we have 

5>*(B fc )*v = £G£)' c^E^'^r^EM^r. (3.18) 



In the last equality we have identified r = — and 
g = — r*. We then see that computing the spectrum of 
fx* by covering its support with boxes of given size is the 
same as computing the spectrum of /i by means of a cov- 
ering of boxes of given mass. That is, one measure is the 
inverse of the other, in the sense of [M. Specializing to 
boxed of fixed size of mass, we can state that computing 
the fixed-mass spectrum of the naive measure fj, on the 
sequence T(N) amounts to the computation of the fixed- 
size spectrum of the dual measure fi* on the approximate 
Cantor set Ct(N), and the other way around. 

In the remaining of this paper we will focus mainly 
on the spectrum of the dual measure fx* for the time 
sequence considered (dual spectrum), as opposite to the 
spectrum of the naive measure (naive spectrum). There- 
fore, in order to alleviate notation we will denote this par- 
ticular dual spectrum and associated magnitudes with- 
out the explicit star-superindex notation, unless other- 
wise stated. 



IV. NUMERICAL RESULTS FOR SYNTHETIC 
TIME SEQUENCES 

In this section we present our estimates for the multi- 
fractal spectrum of some synthetic (computer generated) 
time sequences. First we check our algorithm with two 
measures of known multifractal spectrum. Finally, we 
study the special case of a random signal whose values 
are distributed according to a truncated power law. 

The numerical procedure for computing estimates of 
dimensions D(q) implies the quenched average of the par- 
tition sum over an ensemble of statistically independent 
realizations of the signal, each one with the same length 
N. By quenched averages we refer to the mean value 
of the logarithm of the partition sum, (log Z £ (q)) . As it 
is well-known, this kind of average is more stable and 
less subject to a particular sampling of scarce signifi- 
cance than the annealed average, which would consider 
the logarithm of the mean value of the partition sum, 
log(Z e (g)). In order to obtain results comparable in a 
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straight forward way for any value of q, the linear re- 
gressions to estimate D(q) are always performed over 
the same scaling interval, independently of the particular 
value of q considered. 



A. Uniform random sequence 



Firstly, we analyze a uniform random sequence 
TZ(N,m,a), where the different values t n are uniform 
uncorrelated random variables with mean value m and 
standard deviation a. For our numerical experiments 
we choose m = 100 and a = 10. For a smooth signal 
like 1Z(N, m, a) we expect to obtain a flat multifractal 
spectrum, that is, generalized dimensions equal to unity 
for both naive and dual measures. This expectation is 
confirmed by our computations, which yield generalized 
dimensions satisfying \D{q) — 1| < 0.001 for \q\ < 10, and 
dimensions very close to 1 for 10 < \q\ < 40 



B. Self-similar deterministic sequence 



of len gth N = 10000, together with the analytic spec- 
trum (4.1). The figure shows an excellent agreement be- 
tween our estimates and the expected analytic result, in 
the whole interval of values of q considered, both posi- 
tive and negative. The accuracy of the fit can be slightly 
improved by increasing the sequence length, but the es- 
timates are already quite stable and correct for the value 
of N showed in the figure. Computations performed for 
the naive spectrum yielded an equally good agreement 
with the analytical result. 
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We can construct a fully multifractal sequence starting 
from any self-similar deterministic multifractal measure 
on 1R Jl^,^2|. We considered a fixed-size measure with 
contraction factor r = 1/2 and probabilities p\ = 0.3 
and P2 — 0.7 |32]| and constructed a non- normalized ap- 
proximation of the measure composed by 1.1 x 10 7 points 
by means of the 'Chaos Game' The multifractal se- 
quence was eventually constructed by binning the sample 
points in 5 x 10 4 boxes covering the interval ]0, 1] over 
which the original measure was defined. The value t n 
of the sequence is then given by the occupation number 
of the n-th box. Fig 2(a) depicts such a sequence. Its 
self-similarity seems obvious even to the naked eye. 

The analytical dual spectrum of the sequence is given, 
as a function of the parameter s £ 1R, by the expres- 
siongfH 



log(pf +p s 2 ) 
log(s) 



(4.1) 



D(q s 



1 



log(Pi+P|) 
log(s) 



(Recall that D(q) stands now the fixed-mass spectrum of 
the original naive multifractal measure. The expression 
for its fixed-size spectrum, commonly found in the litera- 
ture, is rather less complex.) In our computations we av- 
eraged over 10 different approximations of the sequence. 
Linear regressions were performed over an interval of 2.5 
decades. Error bars correspond to statistical errors from 
the regression algorithm. In Fig 2(b) we have plotted our 
numerical estimates of the dual spectrum for sequences 
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Fig. 2: a) Succession of 50000 values from a determinis- 
tic multifractal time sequence. Its parameters are r = 1/2, 
pi = 0.3, andp2 = 0.7 (see text), b) Mathematical dual spec- 
trum of the above sequence (full line); the points represent 
our numerical estimates. 



C. Power-law random sequence 

The sequence of transit times seems to be distributed 
according to a truncated power-law of the form 



, . _ /const. t e [0,t [ 
P ^ to >-\at^ te[t ,oo[ 



(4.2) 



(see Eq. (1.1)). In order to explore the applicabil- 
ity of our algorithm to a power-law sequence, we have 
constructed and analyzed a synthetic random sequence 
V(N,to) — {t n } n= i jvi m which each t n is a random 
variable sorted according to the density (|4.2]). In order 
to get results comparable to those of the transit times 
sequences, we will only allow values of x m the range 



G 



2 < x < 3. For the purposes of our computer calcula- 
tions, we generate a synthetic sequence by sampling TV 
values t n according to the rule 



to 



l-1o 



-V(x-i) 



if < i?o 
if 7? > 770 



(4.3) 



where 77 is a uniform random number in ]0, 1] and 770 = 
1 — 1/X- (See Appendix for details.) Given that each term 
t n of any particular realization of the sample depends lin- 
early on /jo, we infer that the multifractal spectrum of the 
sequence will be independent of the particular cut-off to 
chosen. We will report results on Dj^(q), the multifractal 
spectrum computed for an ensemble of sequences of fixed 
length TV. 

When computing the spectrum for any given value of 
X £]2, 3[, we find that for any fixed TV, the results for dif- 
ferent samples of the sequence do not collapse onto the 
same function, but are widely scattered around some av- 
erage position. We explain that effect by the fact that, 
by construction, the signal t n has no upper bound, so 
that it is possible to find that, just perchance, we have 
generated a sample with a particular term t p extremely 
large, in comparison with the expected average maximum 
value (Tm) (that is, a rare event). Is easy to show (see 
Appendix) that in a sequence of TV random variables dis- 
tributed according to a truncated power-law, the average 
maximum value expected scales in the limit of large TV 
as 



(TAf^VV 1 ^- 1 ). 



(4.4) 



In order to get rid of the effect of those rare events, 
we proceed to compute the spectrum of a depleted se- 
quence, in which all the values t n larger than a threshold 
Tm = toN 1 ^^ 1 ^ have been truncated to the value Tm- 
By using this trick, we obtain stable results for all se- 
quence lengths, collapsing onto the same average curve, 
within the error bars. In order to check that our particu- 
lar selection of the threshold does not have an exceedingly 
strong effect on the computed spectra, we have repeated 
our calculations for different values of Tm, finding always 
the same behavior for the generalized dimensions, even 
for a threshold as large as t N. In the computations re- 
ported here, we average for each sequence length over an 
ensemble of 25 different realizations. Linear regressions 
were performed on intervals of 2 decades. Statistical er- 
ror bars are all smaller than 0.01. 

First of all, we observe that for q < 0, the dual spectra 
are always ill-defined, suffering from unacceptable corre- 
lation coefficients and therefore being meaningless. This 
fact seems to be very natural, since, as it is well-known, 
fixed-size algorithms render bad results for negative q. 
However, recall that what we are actually measuring is 
the fixed-m ass spectrum of the naive measure defined in 
Sec. Ill A, so that the fixed-size spectrum of that very 
measure turns out to be well behaved for negative q, and 
ill-defined for positive q, against all previous intuition. 



The reason of this fact is the following: For negative q 
the partition function is dominated by the sparse regions 
of the measure and, for positive q, for the dense regions. 
The bad behavior for q < is a reflection of the presence 
on holes in the support of the dual measure, the only 
source of boxes with abnormally small measure. Going 
back to the naive measure, this means that this measure 
is dominated by a background of a few points with an 
extremely large measure (corresponding to the holes in 
the dual measure), which cause the break down of the al- 
gorithm for positive q. We claim therefore that the dual 
measure as defined in Sec. [II B is the most appropriate 
to characterize extremely non-homogeneous series, like 
the power-law distribution under consideration. 
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Fig. 3: Multifractal dual spectrum for a power-law time 
sequence with exponent \ = 2.22. From top to bottom, 
the curves depict the spectrum for sequences of length 10 6 , 
3 x 10 5 , 10 s , 3 x 10 4 , 10 4 , 3 x 10 3 , and 10 3 , respectively. 

In the range q > 0, for every value of x analyzed we ob- 
serve stable dual spectra, dependent on TV, for TV > 1000. 
When increasing the value of TV, the spectrum becomes 
flatter and flatter. That is to say, the "multifractal- 
ity" of the sequence becomes smaller and smaller, with 
Dn{i) — > 1 for any q, when TV — > 00. This result is 
shown in Fig. 3. A measure of the degree of multifrac- 
tality (multifractality strength), of the sequence could be 
the expression 1 — 0^(00), where 



D N (oo) 



lim D N (q). 

q— >+oo 



(4.5) 



We have computed D^(oo) from linear regressions of the 
partition function computed for a value of q large enough 
to ensure the stability of the estimates. Numerically we 
find that the multifractality strength is a power-law func- 
tion of TV, with an exponent dependent on x 



l-Djv(oo) 



TV 



-7(x) 



(4.6) 



In Fig. 4 we have plotted 1 — Dn(oo) versus TV in log- log 
scale, for different values of x- The change in the slope 
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is evident. In Fig. 5 we represent the estimated values of 
7 as a function of x- It is very well approximated by a 
linear relationship "f(x) ~ X- Our numerical estimates of 
the coefficients of this relation are 



7(x) = (0.48 ± 0.01)x - (0.82 ± 0.02) 



(4.7) 



In the limit of infinite N we will find a flat spectrum 
(uniform measure); however, for any finite value of N 
the deterministic sequences are fully multifractal. 
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Fig. 4: Plot of 1 — Djv(oo) as a function of N for 9 values of \\ 
from top to bottom, \ varies from 2.1 to 2.9, in steps of 0.1. 
The full lines are linear fittings to the power-law behavior. 
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Eq.(4.6) suggests the possibility of some sort of finite- 
size scaling for the multifractal spectrum: We can rewrite 
(4.6) in the form 



l-Av(oo) 
AT--y(x) 



const., 



(4.8) 



that is, in the limit q — > oo, the spectra scales as a power 
law of the sequence length. In view of this last formula, 
one would be tempted to extend the scaling to all val- 
ues of q, defining a renormalized spectrum through the 
expression 



1 - D N (q) 

AT-7(x) 



l-D R (q). 



(4.9) 



The renormalized spectrum Dn(q) is a universal func- 
tion, independent of the length N. It is an intrinsic 
property of the initial time sequence, independent of any 
particular sample, and it can be therefore regarded as its 
true spectrum. 
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Fig. 6: Finite-size scaling of the multifractal dual spectrum 
for power-law time sequence with exponent \ = 2.22 



In Fig. 6 we have tested the scaling ansatz ( |4.9| ) for 
actual computations. The best collapse is achieved for 
sequences with length in between 10 4 and 10 6 , and for 
an exponent 7' = 0.265. The power-law sequence consid- 
ered has a distribution exponent x — 2-22 and a predicted 



value 7 = 0.25 according to Eq. (4.7), quite close to the 
actual value. 



NUMERICAL RESULTS FOR TRANSIT 
TIMES SEQUENCES 



We now turn to the numerical analysis of the sequence 
of SOC transit times S(N, L). By construction, the value 
T n is the time spent into the pile by the n-th grain in a 
series of N consecutive throws. It is conceivable that 
the landing of a tracer may provoke an avalanche which 
would eventually evacuate out of the pile the very tracer 
that caused it. In such a case, we assign a value T = 1 
to the transit time of that particular tracer. We have 
therefore T n e [l,oo[. Since the computer time devoted 
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to any simulation is always a finite amount, one has to 
stop the run at some point, leaving inside the pile, with 
nonzero probability, some of the tracers thrown at inter- 
mediate stages of the simulation. These tracers which 
did not emerge at the end of the run would represent a 
gap in the sequence S(N, L). We fill these gaps by shift- 
ing the sequence one site to the left at the points n when 
a tracer did not come out. We have also considered se- 
quences in which each gap was filled with a lower bound 
of its corresponding transit time, estimated by substract- 
ing the time of addition of tha gap to the total time that 
the simulation was running. The results obtained with 
both procedures were identical, within the error bars. 
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Fig. 7: Multifractal dual spectrum for SOC sequences from 
a ricepile of size L = 100. The different plots correspond 
to different sequence lengths; from top to bottom, N = 10 5 , 
3 x 10 4 , 10 4 , 3 x 10 3 , and 10 3 . 
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Fig. 8: Multifractal dual spectrum for SOC sequences of 
length N = 10 4 . The different plots correspond to different 
system sizes; from top to bottom, L — 25, 50, 100, 200, 400, 
800, and 1600. 

We work with sequences of total length M — 10 6 points 



from simulations of the one dimensional Oslo model of 
size L = 25, 50, 100, 200, 400, 800, and 1600. In order to 
average our partition sum, we proceed to decompose the 
sequences into subsequences of length N <C M, and per- 
form the averages over the sample of the resulting M/N 
subsequences. When computing the spectra, however, we 
find that they do not stabilize well. This is again due to 
the presence of rare events: In a subsequence of length N 
there are some points with extremely large relative mea- 
sure, corresponding to tracers which spent a long time 
inside the pile. In order to correct this effect, we proceed 
in the same way as we did in the random power-law signal 
above: We truncate the larg est e vents up to a maximum 
cut-off T m ■ In view of Eqs. ([□]) and the SOC sig- 

nal is akin to a truncated power- law distribu ted sequence 
with cut-off to ~ L"; comparing with Eq. (4.4), we se- 
lect T M = L v N l /&- x \ with v = 1.30 and X = 2.22, 
according to the simulations. Our results are the spectra 
Z3jv.l(<?)j computed for an ensemble of sequences of fixed 
length N, coming from a ricepile of size L. 

With the expertise we gained from the analysis of the 
random power-law signal, we would expect the multi- 
fractal spectrum of any SOC sequence to be ill-defined 
for q < 0, to depend on the length N, and to be indepen- 
dent of the cut-off, that is, of the system size L. The first 
prediction turns out to be true; for q < the poor correla- 
tion coefficients yield meaningless estimations. However, 
for q > we obtain stable spectra depending on both N 
and L. They show an even more striking property; the 
spectra decrease monotonically (become flatter) with N 
and increase (become steeper) with L. This behavior is 
shown in Figs. 7 and 8. 
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Fig. 9: Finite-size scaling of the multifractal dual spectrum 
for SOC sequences. 

In a similar way as we did for the synthetic signal, we 
proceed to investigate the degree of multifractaly of the 
SOC sequence. Studying the same strength parameter, 
we find that the magnitude 1 — Dn,l{oo) can be fitted 
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as a double power-law, both in N and L, that is, 
1 - D N , L (oo) ~ N-^L^ 2 . 



(5.1) 



Our estimates are 71 = 0.27 ± 0.02 and 72 = 0.32 ± 0.02. 
These results are valid in the range N > 10000 and 
L < 400. 

The previous formula suggests again the possibility of 
constructing a renormalized spectrum, universal for all 
values of q, and independent of N and L. This is done 
by plotting the finite-size relationship 



l-D NtL {q) 



l-D R (q). 



(5.2) 



The validity of this scaling is checked up in Fig. 9. The 
plotted spectra correspond to the smaller values of L and 
larger values of N for which the relation (5.1) holds. The 
best collapse is obtained for effective exponents j[ = 0.29 
and 7 2 = 0.34, very close to the ones predicted in the 
limit q — > 00. The rescaled spectra collapse onto a unique 
function, which is interpreted again as a renormalized 
spectrum, in the sense that it is a property of the intrin- 
sic dynamics of the ricepile where the data came from, 
and independent of particular samples considered when 
computing it. 

This scaling behavior can be accounted for by the effect 
of the correlations inside the SOC sequence. No depen- 
dence whatsoever on the system size (the cut-off) was 
observed in the synthetic power-law distributed signal in 
Sec. [VC. The only difference between that signal and 
the SOC one resides in the correlations. While the dif- 
ferent points in the synthetic sequence are completely 
uncorrelated by construction, the SOC transit times suf- 
fer obviously from long-range correlations. This fact is 
easy to realize when one considers that grains introduced 
into the pile at widely scattered initial times can emerge 
at the same instant in a single gigantic avalanche. 

As a numerical experiment, we have estimated the cor- 
relation length in our SOC sequences as the minimum 
length N above which an R/S analysis |53 provides a 
Hurst exponent close to 0.5. Our estimates show that for 
L < 400 the sequences become roughly uncorrelated for 
lengths larger that N = 10 4 , whereas no serious estimate 
can be done for L > 400. This result seems to be in 
contradiction with our multifractal scaling, since in the 
range of validity of Eq. fl5.l|), the R/S analysis predicts a 
complete decorrelation and, hence, an independence on 
the system size. We interpret our results as a hint to- 
wards the existence of more deep intrinsic correlations 
that those revealed by a simple R/S analysis. 



VI. CONCLUSIONS 

In this paper we have investigated the multifractal 
properties of sequences of transit times of individual 
grains inside an Oslo ricepile. To this purpose, we have 
developed a fixed-mass multifractal algorithm, yielding 



the so-called dual spectrum, particularly well-suited to 
deal with highly inhomogeneous one dimensional mea- 
sures (in our case, time series). This is particularly for the 
transit time sequences, which are power-law distributed 
and are hence constituted at any length scale by a more 
or less average flat background, interspersed by relatively 
infrequent huge peaks. 

The main result of our analysis is the finite-size scaling 
relation (5.2). This scaling shows a particular behavior: 
The dual spectrum tends to decrease when increasing the 
sequence length N, whereas it tends to increase with the 
systems size L. While the first statement is in complete 
agreement with numerical experiments on synthetic un- 
correlated power-law sequences, the second constit utes a 
completely unexpected result: As we show in Sec. IV C , 
the spectra of an uncorrelated random power law signal 
do not depend on the distribution's cut-off. Since the cut- 
off is related to the system size of the ricepile, we should 
expect in the SOC case to obtain results independent of 
L. That is not the case, however, in our computations. 

The renormalized spectrum defined in (|5.2| ) allows one 
to get rid of those finite-size effects, and constitutes a 
magnitude that can be associated to the very ricepile dy- 
namics, not influenced by the hazards of the samples used 
in its estimation. 

We interpret the initial L dependence as an effect of 
the extremely long correlations in the transit time se- 
quence. As the authors point out in Ref. pQ], the fact 
that the average speed of the tracers decreases with the 
system size proves that there are correlations all along 
the system. These correlations show up even more spec- 
tacularly when analyzing the multifractal properties of 
the sequences. A simple R/S analysis seems to show an 
absence of correlations for L < 400 and N > 10 4 . Hence, 
it could seem reasonable that, for these values of the pa- 
rameters, the spectra should become independent of L. 
This is not the case, however. We conclude, therefore, 
that the transit time sequences indeed posess correlations 
of a range far larger than that possibly revealed by the 
R/S analysis, correlations which are made evident only 
in our more sophisticated multifractal analysis. 
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APPENDIX: 

In this appendix we derive some useful properties of a 
truncated power-law random variable. Consider a ran- 



dom variable t distributed according to the density (4.2) 



Continuity of the density at t = to imposes the actual 
form 



a t n 



te [o,to[ 

t e [to,oo[ 



(Al) 



If x > 1) then the density is normalizable, with a nor- 
malization constant 



t dt 
to to 



dt 



X 



X 



1 



(A2) 



If we demand that the first moment of the distribution 
does exist, we have to impose x > 2 to obtain 

(t) = J t p(t, t ) dt = —t ~ t , (A3) 

where A = (x - 2)/(x - 1). 

The distribution function P(t,to) = J p(t,ta)dt has 
the form 



P(t,t ) 



x-i t 

X t 



1 



X (to) 



-x+i 



is [0,t [ 
t e [to)°o[ 



(A4) 



In order to sample a sequence according to this distri- 
bution, we use the inversion method ]3q |: We equate 
the distribution function to a uniform random number 
r\ and obtain the corresponding value of t by inverting 
Pit, to) = n. It i s ea sy to check that the resulting sample 
is given by Eq. (4.3). 

Consider now that we sort N independent random vari- 
ables according to the distribution p(t,to), obtaining the 
sample {ti, . . . tjv}- Define Tm as the maximum value in 
this particular sample, Tm = max{ti, . . .iAr}. We want 
to comp ute the average value (Tm), weighted with the 
density ( Al). It is easy to see that the probability of this 
maximum value being lesser or equal than Tm is just 
equal to the probability of all the individual values t n 
being on their turn lesser or equal than Tm ■ This means 
that the distribution function of the maximum value Tm 
is just 



U{T M ,N) = P{T M ,t ) 



N 



(A5) 



By differentiating Eq. ( ]A5[ ) we get the probability density 
of maximum values 



n(T M ,N) 



dH(T M ,N) 



dT, 



A I 



N 
N 



m (* 



N-l 



-x+i 



N-l 



u i 



T M e [0,* [ 



T M 6 [to, oof 



(A6) 



The average maximum value that we expect to observe in N samples of the initial power law distribution will then be 



(Tm) 



Tm tt(Tm , t ) dT M - 



After substituting Eq. ( |A6| ), we obtain 

(T M ) N 



to 



N-l 



X 



N 



N 
Jo 



i - -e /x 

X 



N-l 



dt 



(A7) 



(A8) 



In the limit N — > oo, the only contribution in the last integral comes from values of £ very close to 0. We can therefore 
evaluate the leading behavior for large N by expanding the integrand in Taylor series, keeping only the first order: 



x 



N-l 



d£ = J exp /(iV-l)ln^l- ~£ 
~ J exp|-(JV-l)^J^ 



N-l 

X 



r(A). 



(A9) 



In estimating the last integral we have extended to infinity the upper limit, approximation allowed again in the limit 
of large N. 

Collecting everything we get finally 



(Tm) 
to 



N 
N - 1 



exp < —N In 



X 



+ x r(A)JV(JV - I) - 



(A10) 
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The first term decays exponentially. Hence in the limit of large N, the leading behavior is given by 

(T M )~t N 1 - x = t N 1 /(x- 1 \ (All) 
up to a constant prefactor, depending only on \. 
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